Raster vs Vector
Install packages we will need.
Module 3 for William Norfolk
library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.0 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.3 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:tidyr':
##
## extractlibrary(sp)
library(mapview)# install.packages("ggplot2")
# install.packages("sp")
# install.packages("raster")
# install.packages("rgdal")
# install.packages("brew")
# install.packages("systemfonts")
# install.packages("mapview")
# install.packages("spatstat")
# install.packages("tmaptools")
# install.packages("elevatr")
# install.packages("rgbif")
# install.packages("remotes")
# remotes::install_github("Nowosad/spDataLarge")
# install.packages("spData")
# install.packages("gstat")
# install.packages("fields")
# install.packages("foreign")
# install.packages("velox")As we’ve seen, vector data at it’s simplest is a set of points with attributes, and can be stored as easily as a table as it can a shapefile or other spatial data format.
library(raster)
dbasin_pts<-shapefile("../IDEAS_Spatial_Ecology_Workshop/dbasin_pts.shp")We made the shapefile above from a .csv that we downloaded. We can convert it back to a .csv without any loss of information.
write.csv(dbasin_pts@data,"./dbasin_pts.csv") #Saved in local spatial folderVector lines are made up of points and connecting arcs. Lines can have values that indicate their width (number of lanes for roads, bank to bank for streams) and can have fields indicating direction – allowing for topology and network analysis (e.g. if bridge X is destroyed, what’s the shortest route from A to B).
Let’s load the stream network for Athens-Clarke county from the USGS National Hydrography Dataset.
library(spData)
library(spDataLarge)
#USGS National Hydrogeography Dataset
download.file(url = "https://www2.census.gov/geo/tiger/TIGER2019/LINEARWATER/tl_2019_13059_linearwater.zip", destfile = "tl_2019_13059_linearwater.zip")
unzip(zipfile = "tl_2019_13059_linearwater.zip")#Another handy function for other types of analysis to note. Command unzips the designated folder and saves contents at the local path
library(raster) # despite the name, includes the basic vector commands
athens_streams<-shapefile("tl_2019_13059_linearwater.shp")
#Note the shapefile command can both load and save shapefiles at the designated pathLet’s view the lines and query them using the mapview package.
#Athens streams
mapview(athens_streams)#Mapview requires spatial data points/lines data frame to produce visuals. Must have .prj file to designate CRSNOTICE that the LINEARID is a unique id for stream reaches that is ordered to indicate direction of flow. (Explanations of attributes are given in the metadata.)
NOTE that the data structure for lines is more complicated than for points. All segments have a FROM point, a TO point and a length. Lines are made up of a series of connected nodes. In vector line, spatial resolution is in part a function of number of nodes. More nodes more resolution.
We can buffer lines or points (or polygons) to make polygons. To do this, we need to project the data into planar coordinates to get planar distance units. We’ll use UTM, but zone 16 this time, since Athens is in eastern Georgia.
library(sp)
#adds distance units to the coordinates to form polygons that define a designated area
#Below we transform the data by changing the CRS designation to zone 16
athens_streams_utm<-spTransform(athens_streams,crs("+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs"))
# Buffering by 10 meters (the units of the projected layer) on either side
# This assigns width
library(rgeos)## rgeos version: 0.5-3, (SVN revision 634)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-1
## Polygon checking: TRUEathens_stream_buff<-gBuffer(athens_streams_utm,byid=TRUE,width=10)
# Because we chose to buffer by id, we end up with as many buffer polygons as stream segments
library(mapview)
mapview(athens_streams_utm) + mapview(athens_stream_buff, col.regions = "red") #Added color scheme to show layers#This adds width to the mapview so now if zoomed in will see more than just a delineated line. Each stream segment (area between two nodes) is now a polygon (sort of elongated oval) and stretches the lenght of the stream.We can also make polygons from bounding boxes.
#Function to create an exten object in R must feed it a raster or matrix
e<-extent(athens_stream_buff)
#Use as() function to designate the object as class SpatialPolygons
athens_box <- as(e, "SpatialPolygons")
#Assign the coordinate system
proj4string(athens_box) <- "+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs"
#Add layer
mapview(athens_box, col.regions = "green") + mapview(athens_streams_utm) + mapview(athens_stream_buff, col.regions = "red") #Adding the third layer back#Renders a polygon over the enture basin supplied in the athens_stream_bluff dataathens_box has no attributes. Simple to add them. Let’s make a column called “Name” with the value “Athens bounding box”.
athens_box$Name<-"Athens bounding box"NOTE that polygons are composed of a set of points (nodes) that define segments of the polygon (the more nodes, the higher the spatial resolution), which also has an area and a perimeter. ONLY when polygons are projected into planar coordinates, can we calculate area.
#calculating area
athens_box$km2<-gArea(athens_box)/1000000 # divide by 10^6 to convert m2 to km2Coming full circle, polygons can be reduced to centroids, the point that marks the rough geographic center of the shape. The
library(geosphere) #need geosphere here. Installed to R server
# Get UTM 16 coordinates of athens_stream_buff polygon centroids
head(geosphere::centroid(athens_stream_buff))## [,1] [,2]
## 0 840008.5 3759775
## 1 841574.6 3757348
## 2 843897.1 3756029
## 3 842960.2 3758165
## 4 843005.5 3758796
## 5 833821.5 3764742#Basically assign x and y locations for each polygon. Called the centroidWhile polygons have a complex data structure, it’s efficient and easy to understand. Polygons can also be overlapping as in shapefiles of species ranges, and can be nested within larger regions. AND, all vector data, whether points, lines or polygons, can be linked to endless attributes that are specific to a location.
However, polygons (and vector data in general) are not great at representing continuous values. For this we generally prefer the raster format.
Better to use polygons/vector for discrete items or things that are sparse.
The utility of the raster format may be best appreciated when interpolating data from a set of locations to create a layer of values across a larger area. Here we have data on precipitation in 2019 from weather stations throughout Georgia (source: http://www.georgiaweather.net/).
#Prepping data to visualize spatial distribution of percipitation in an area
# Read in precip data
ga_precip_2019<-read.csv("../IDEAS_Spatial_Ecology_Workshop/ga_precip_2019.csv")
head(ga_precip_2019)## Station.Name X2019 Normal1971_2000 Departure
## 1 Alapaha Berrien Georgia Georgia 36.98 49.72 -12.74
## 2 Alma Georgia 44.56 49.78 -5.22
## 3 Alpharetta Georgia 53.62 46.89 6.72
## 4 Arabi Georgia 64.23 52.04 12.19
## 5 Arlington Georgia 51.34 60.13 -8.79
## 6 Attapulgus Georgia 46.88 50.07 -3.19# Rename column to show that it records precipitation (inches)
colnames(ga_precip_2019)[2]<-"precip2019"The table includes place names (locations of weather stations), average max temp values and coordinates. Let’s get coordinates for these locations.
library(tmaptools) #Package for the geocode family of functions
#Need to change Station.Name to character first
ga_precip_2019$Station.Name <- as.character(as.factor(ga_precip_2019$Station.Name))
#Pull coordinates using geocode_OSM. This function will accept character strings or direct character input so must recode classes to char if using a dataframe
station_coords<- geocode_OSM(ga_precip_2019$Station.Name)# How to do this? geocode_OSM!
head(station_coords)## query lat lon lat_min lat_max
## 1 Alapaha Berrien Georgia Georgia 31.37862 -83.22247 31.37381 31.38694
## 2 Alma Georgia 41.69918 42.48374 41.67918 41.71918
## 3 Alpharetta Georgia 34.07096 -84.27473 34.02821 34.11735
## 4 Arabi Georgia 31.83157 -83.73795 31.80485 31.85153
## 5 Arlington Georgia 31.43990 -84.72492 31.42070 31.45287
## 6 Attapulgus Georgia 30.74908 -84.48380 30.74220 30.75680
## lon_min lon_max
## 1 -83.22705 -83.21664
## 2 42.46374 42.50374
## 3 -84.35944 -84.20122
## 4 -83.75856 -83.70387
## 5 -84.74235 -84.70828
## 6 -84.49236 -84.47538locations<-station_coords[,c(3,2)]Now we turn these into a SPDFe as before, because we need planar coordinates in order to interpolate.
#create a new spatial points dataframe as to designate a new crs to get planar coordinates
precip_pts<-SpatialPointsDataFrame(locations, #needed to determine locations to add to SPDF
ga_precip_2019,
coords.nrs = numeric(0),
proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))Let’s project to UTM 17 which is the preferred projection for statewide data in Georgia.
#project correct crs
precip_pts_utm<-spTransform(precip_pts,crs("+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"))
# write out as shapefile
shapefile(precip_pts_utm,"ga_precip_pts_utm_2019.shp", overwrite = TRUE) #adding overwrite to fix knitting issue## Warning in rgdal::writeOGR(x, filename, layer, driver = "ESRI Shapefile", :
## Field names abbreviated for ESRI Shapefile driverA raster has a coordinate system, and a cell-size (resolution) that along with its extent, determines the number of columns and rows of cells. Let’s use the extent (bounding box) of another existing UTM 17 layer (Conservation Lands) as a template.
#extent = boundry limits of a defined area we used this to make the Athens box above
# can pull the extent directly from the shapefile usinf raster::extent
# Bring in the conservation lands
conlands<-shapefile("../IDEAS_Spatial_Ecology_Workshop/ConservationLands_2019.shp")## Warning in .local(x, ...): .prj file is missing# Get conlands extent
e<-raster::extent(conlands)Now let’s make an empty raster with 1 km2 cells over the extent of the Conservation Lands.
#res = number of cells which in turn determines the "resolution" of the raster
precipr<-raster(ext=e,res=1000) #the plain raster command makes a new raster set by the bounds of the extent file "e" and divided into pixels by the resolution value
#need to add a .proj
proj4string(precipr)<-"+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"
precipr## class : RasterLayer
## dimensions : 500, 431, 215500 (nrow, ncol, ncell)
## resolution : 1000, 1000 (x, y)
## extent : 83446.45, 514446.4, 3379976, 3879976 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0NOTICE that this blank raster has 500 rows and 431 columns. Now let’s interpolate precipitation (inches) in 2019, average annual precipitation from 1971 to 2000, and departure from that average in 2019 over this grid of cells using a thin plate spline model.
Are the number of rows and columns determined by the shape of the extent?
library(fields) #additional package for curves and other spatial tools## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.5-1 (2019-12-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: maps##
## Attaching package: 'maps'## The following object is masked from 'package:purrr':
##
## map## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code#Tps function: Fits a thin plate spline surface to irregularly spaced data
precip <- Tps(coordinates(precip_pts_utm), precip_pts_utm$precip2019)
#makes a raster layer with interpolated values using a model object
precip2019_tps <- interpolate(precipr, precip)
normal<- Tps(coordinates(precip_pts_utm), precip_pts_utm$Normal1971_2000)
normal1971_2000_tps <- interpolate(precipr, normal)
# view the results
mapview(precip2019_tps) # Notice that we can query cells in mapviewmapview(normal1971_2000_tps) # Notice that we can query cells in mapview#Big differene in the percipThis is somewhat coarse, but shows the general trend. Spatial statistical approaches are beyond the scope of this workshop. The purpose here is simply to demonstrate the raster model and its usefulness for mapping continuous values. Let’s look at the layers.
precip2019_tps## class : RasterLayer
## dimensions : 500, 431, 215500 (nrow, ncol, ncell)
## resolution : 1000, 1000 (x, y)
## extent : 83446.45, 514446.4, 3379976, 3879976 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : layer
## values : 39.00672, 80.55721 (min, max)normal1971_2000_tps## class : RasterLayer
## dimensions : 500, 431, 215500 (nrow, ncol, ncell)
## resolution : 1000, 1000 (x, y)
## extent : 83446.45, 514446.4, 3379976, 3879976 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : layer
## values : 40.9473, 85.47134 (min, max)#to layer rasters do all need to have the same demension, resolution, extent, and crs? I would imagine so...if data is missing from one layer can it be added/interpolates/labeled as NA?Notice that this raster has a range of values for precipitation (~39-81 inches), but no other attributes. We can link labels to the precip values, but not to cell locations as in vector GIS. But, while the data model is simpler, the files (especially when high resolution) are generally much larger than vector files of the same extent, because they include values for, in this case, 215500 cells. If we had more weather stations, it might have made sense to interpolate to a raster with even finer resolution.
Raster is limited by computing power at high resolution.
We can create a histogram of 2019 precip values for Georgia.
#cannot use ggplot unless data is converted to dataframe
hist(precip2019_tps)
NOTE that raster data often takes up a LOT of memory and therefore is not usually stored in .RData files. So if you don’t want to redo the steps, you need to write out the raster.
# write to file as a GeoTiff
#writeRaster(precip2019_tps,"precip_ga_2019_tps.tif") #write a Geotiff into local path
# to read back in
#precip2019_tps<-raster("precip_ga_2019_tps.tif")
#will leave this code chunk commented out since the file is written now and will be created baove when knittingSee if you can repeat these steps for the drought (Departure) column of the precip2019 Georgia data.
#First, lets make a new raster
departr<-raster(ext=e,res=1000)#will keep the same extent = e
#need to add a .proj
proj4string(departr)<-"+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"#next we pul Departure from precip_pts_utm and interpolate as above
depart <- Tps(coordinates(precip_pts_utm), precip_pts_utm$Departure)
depart_tps <- interpolate(departr, depart)
# view the results
mapview(depart_tps) #raster layer of droughts#view histogram
hist(depart_tps)
#Lastly we can write the file out to the local path and pull it back in if needed
#writeRaster(depart_tps,"depart_tps.tif") #write a Geotiff into local path
# to read back in
#depart_tps<-raster("depart_tps.tif")
#will leave this code chunk commented out since the file is written now and will be created baove when knittingWe’ve just looked at a continuous raster. Not all continuous rasters derive from interpolated values. In many cases, we have sensors that provide continuous coverage of for example, night-time lights. But, landcover classifications (from remotely sensed data) are usually served as raster data. In thise case, the data is discrete or categorical with numbers standing in for particular landcover classes.
Think continious monitoring data from a weather station, need to “bin” the values to define a specific raster category as a range of these values and assign the appropriate color/demarcation/etc. for the raster to display
Let’s look at the most recent (2015) land use classification of Georgia done by the Georgia Land Use Trends project (http://narsal.uga.edu/glut/) and available from the Georgia GIS Data Clearinghouse (https://data.georgiaspatial.org/).
ga_lc<-raster("../IDEAS_Spatial_Ecology_Workshop/glut_lc_2015.img")
ga_lc## class : RasterLayer
## dimensions : 18910, 18791, 355337810 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 7515, 571245, 3342615, 3909915 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : /opt/home/wn45960/IDEAS_Spatial_Ecology_Workshop/glut_lc_2015.img
## names : glut_lc_2015
## values : 0, 93 (min, max)In this case, the raster is accompanied by a .dbf that gives us the a description of the classes.
library(foreign) #Package used to read in and write out data in atypical formats
# read database file
vat<-read.dbf('../IDEAS_Spatial_Ecology_Workshop/glut_lc_2015.img.vat.dbf')
vat## Value Count Red Green Blue Opacity Class
## 1 0 184160185 0 0 0 255 Background
## 2 1 0 0 0 0 255 <NA>
## 3 2 0 0 0 0 255 <NA>
## 4 3 0 0 0 0 255 <NA>
## 5 4 0 0 0 0 255 <NA>
## 6 5 0 0 0 0 255 <NA>
## 7 6 0 0 0 0 255 <NA>
## 8 7 185642 249 249 219 255 Beach/Dune/Mud
## 9 8 0 0 0 0 255 <NA>
## 10 9 0 0 0 0 255 <NA>
## 11 10 0 0 0 0 255 <NA>
## 12 11 4216554 91 91 209 255 Open Water
## 13 12 0 0 0 0 255 <NA>
## 14 13 0 0 0 0 255 <NA>
## 15 14 0 0 0 0 255 <NA>
## 16 15 0 0 0 0 255 <NA>
## 17 16 0 0 0 0 255 <NA>
## 18 17 0 0 0 0 255 <NA>
## 19 18 0 91 91 91 255 <NA>
## 20 19 0 0 0 0 255 <NA>
## 21 20 0 107 142 162 255 <NA>
## 22 21 0 0 0 0 255 <NA>
## 23 22 15115798 249 137 158 255 Low Intensity Urban
## 24 23 0 0 0 0 255 <NA>
## 25 24 3807234 209 68 82 255 High Intensity Urban
## 26 25 0 0 0 0 255 <NA>
## 27 26 0 0 0 0 255 <NA>
## 28 27 0 0 0 0 255 <NA>
## 29 28 0 0 0 0 255 <NA>
## 30 29 0 0 0 0 255 <NA>
## 31 30 0 0 0 0 255 <NA>
## 32 31 6046366 249 249 130 255 Clearcut/Sparse
## 33 32 0 0 0 0 255 <NA>
## 34 33 0 249 160 0 255 <NA>
## 35 34 195116 185 185 185 255 Quarries/Strip Mines/Rock Outcrop
## 36 35 0 0 0 0 255 <NA>
## 37 36 0 0 0 0 255 <NA>
## 38 37 0 0 0 0 255 <NA>
## 39 38 0 0 0 0 255 <NA>
## 40 39 0 0 0 0 255 <NA>
## 41 40 0 0 0 0 255 <NA>
## 42 41 35023511 175 223 153 255 Deciduous Forest
## 43 42 38301106 58 153 86 255 Evergreen Forest
## 44 43 6168904 108 193 108 255 Mixed Forest
## 45 44 0 0 0 0 255 <NA>
## 46 45 0 0 0 0 255 <NA>
## 47 46 0 0 0 0 255 <NA>
## 48 47 0 0 0 0 255 <NA>
## 49 48 0 0 0 0 255 <NA>
## 50 49 0 0 0 0 255 <NA>
## 51 50 0 0 0 0 255 <NA>
## 52 51 0 0 0 0 255 <NA>
## 53 52 0 0 0 0 255 <NA>
## 54 53 0 0 0 0 255 <NA>
## 55 54 0 0 0 0 255 <NA>
## 56 55 0 0 0 0 255 <NA>
## 57 56 0 0 0 0 255 <NA>
## 58 57 0 0 0 0 255 <NA>
## 59 58 0 0 0 0 255 <NA>
## 60 59 0 0 0 0 255 <NA>
## 61 60 0 0 0 0 255 <NA>
## 62 61 0 0 0 0 255 <NA>
## 63 62 0 0 0 0 255 <NA>
## 64 63 0 0 0 0 255 <NA>
## 65 64 0 0 0 0 255 <NA>
## 66 65 0 0 0 0 255 <NA>
## 67 66 0 0 0 0 255 <NA>
## 68 67 0 0 0 0 255 <NA>
## 69 68 0 0 0 0 255 <NA>
## 70 69 0 0 0 0 255 <NA>
## 71 70 0 0 0 0 255 <NA>
## 72 71 0 0 0 0 255 <NA>
## 73 72 0 0 0 0 255 <NA>
## 74 73 0 210 22 210 255 <NA>
## 75 74 0 0 0 0 255 <NA>
## 76 75 0 0 0 0 255 <NA>
## 77 76 0 0 0 0 255 <NA>
## 78 77 0 0 0 0 255 <NA>
## 79 78 0 0 0 0 255 <NA>
## 80 79 0 0 0 0 255 <NA>
## 81 80 0 197 146 104 255 <NA>
## 82 81 35830800 211 160 108 255 Row Crop/Pasture
## 83 82 0 0 0 0 255 <NA>
## 84 83 0 223 170 119 255 <NA>
## 85 84 0 0 0 0 255 <NA>
## 86 85 0 0 0 0 255 <NA>
## 87 86 0 0 0 0 255 <NA>
## 88 87 0 0 0 0 255 <NA>
## 89 88 0 0 0 0 255 <NA>
## 90 89 0 0 0 0 255 <NA>
## 91 90 0 0 0 0 255 <NA>
## 92 91 24332457 172 73 137 255 Forested Wetland
## 93 92 1531586 193 160 209 255 Non-Forested Salt/Brackish Wetland
## 94 93 422551 142 130 223 255 Non-Forested Freshwater Wetland
## 95 94 0 0 0 0 0 <NA>
## 96 95 0 0 0 0 0 <NA>
## 97 96 0 0 0 0 0 <NA>
## 98 97 0 0 0 0 0 <NA>
## 99 98 0 0 0 0 0 <NA>
## 100 99 0 0 0 0 0 <NA>
## 101 100 0 0 0 0 0 <NA>
## 102 101 0 0 0 0 0 <NA>
## 103 102 0 0 0 0 0 <NA>
## 104 103 0 0 0 0 0 <NA>
## 105 104 0 0 0 0 0 <NA>
## 106 105 0 0 0 0 0 <NA>
## 107 106 0 0 0 0 0 <NA>
## 108 107 0 0 0 0 0 <NA>
## 109 108 0 0 0 0 0 <NA>
## 110 109 0 0 0 0 0 <NA>
## 111 110 0 0 0 0 0 <NA>
## 112 111 0 0 0 0 0 <NA>
## 113 112 0 0 0 0 0 <NA>
## 114 113 0 0 0 0 0 <NA>
## 115 114 0 0 0 0 0 <NA>
## 116 115 0 0 0 0 0 <NA>
## 117 116 0 0 0 0 0 <NA>
## 118 117 0 0 0 0 0 <NA>
## 119 118 0 0 0 0 0 <NA>
## 120 119 0 0 0 0 0 <NA>
## 121 120 0 0 0 0 0 <NA>
## 122 121 0 0 0 0 0 <NA>
## 123 122 0 0 0 0 0 <NA>
## 124 123 0 0 0 0 0 <NA>
## 125 124 0 0 0 0 0 <NA>
## 126 125 0 0 0 0 0 <NA>
## 127 126 0 0 0 0 0 <NA>
## 128 127 0 0 0 0 0 <NA>
## 129 128 0 0 0 0 0 <NA>
## 130 129 0 0 0 0 0 <NA>
## 131 130 0 0 0 0 0 <NA>
## 132 131 0 0 0 0 0 <NA>
## 133 132 0 0 0 0 0 <NA>
## 134 133 0 0 0 0 0 <NA>
## 135 134 0 0 0 0 0 <NA>
## 136 135 0 0 0 0 0 <NA>
## 137 136 0 0 0 0 0 <NA>
## 138 137 0 0 0 0 0 <NA>
## 139 138 0 0 0 0 0 <NA>
## 140 139 0 0 0 0 0 <NA>
## 141 140 0 0 0 0 0 <NA>
## 142 141 0 0 0 0 0 <NA>
## 143 142 0 0 0 0 0 <NA>
## 144 143 0 0 0 0 0 <NA>
## 145 144 0 0 0 0 0 <NA>
## 146 145 0 0 0 0 0 <NA>
## 147 146 0 0 0 0 0 <NA>
## 148 147 0 0 0 0 0 <NA>
## 149 148 0 0 0 0 0 <NA>
## 150 149 0 0 0 0 0 <NA>
## 151 150 0 0 0 0 0 <NA>
## 152 151 0 0 0 0 0 <NA>
## 153 152 0 0 0 0 0 <NA>
## 154 153 0 0 0 0 0 <NA>
## 155 154 0 0 0 0 0 <NA>
## 156 155 0 0 0 0 0 <NA>
## 157 156 0 0 0 0 0 <NA>
## 158 157 0 0 0 0 0 <NA>
## 159 158 0 0 0 0 0 <NA>
## 160 159 0 0 0 0 0 <NA>
## 161 160 0 0 0 0 0 <NA>
## 162 161 0 0 0 0 0 <NA>
## 163 162 0 0 0 0 0 <NA>
## 164 163 0 0 0 0 0 <NA>
## 165 164 0 0 0 0 0 <NA>
## 166 165 0 0 0 0 0 <NA>
## 167 166 0 0 0 0 0 <NA>
## 168 167 0 0 0 0 0 <NA>
## 169 168 0 0 0 0 0 <NA>
## 170 169 0 0 0 0 0 <NA>
## 171 170 0 0 0 0 0 <NA>
## 172 171 0 0 0 0 0 <NA>
## 173 172 0 0 0 0 0 <NA>
## 174 173 0 0 0 0 0 <NA>
## 175 174 0 0 0 0 0 <NA>
## 176 175 0 0 0 0 0 <NA>
## 177 176 0 0 0 0 0 <NA>
## 178 177 0 0 0 0 0 <NA>
## 179 178 0 0 0 0 0 <NA>
## 180 179 0 0 0 0 0 <NA>
## 181 180 0 0 0 0 0 <NA>
## 182 181 0 0 0 0 0 <NA>
## 183 182 0 0 0 0 0 <NA>
## 184 183 0 0 0 0 0 <NA>
## 185 184 0 0 0 0 0 <NA>
## 186 185 0 0 0 0 0 <NA>
## 187 186 0 0 0 0 0 <NA>
## 188 187 0 0 0 0 0 <NA>
## 189 188 0 0 0 0 0 <NA>
## 190 189 0 0 0 0 0 <NA>
## 191 190 0 0 0 0 0 <NA>
## 192 191 0 0 0 0 0 <NA>
## 193 192 0 0 0 0 0 <NA>
## 194 193 0 0 0 0 0 <NA>
## 195 194 0 0 0 0 0 <NA>
## 196 195 0 0 0 0 0 <NA>
## 197 196 0 0 0 0 0 <NA>
## 198 197 0 0 0 0 0 <NA>
## 199 198 0 0 0 0 0 <NA>
## 200 199 0 0 0 0 0 <NA>
## 201 200 0 0 0 0 0 <NA>
## 202 201 0 0 0 0 0 <NA>
## 203 202 0 0 0 0 0 <NA>
## 204 203 0 0 0 0 0 <NA>
## 205 204 0 0 0 0 0 <NA>
## 206 205 0 0 0 0 0 <NA>
## 207 206 0 0 0 0 0 <NA>
## 208 207 0 0 0 0 0 <NA>
## 209 208 0 0 0 0 0 <NA>
## 210 209 0 0 0 0 0 <NA>
## 211 210 0 0 0 0 0 <NA>
## 212 211 0 0 0 0 0 <NA>
## 213 212 0 0 0 0 0 <NA>
## 214 213 0 0 0 0 0 <NA>
## 215 214 0 0 0 0 0 <NA>
## 216 215 0 0 0 0 0 <NA>
## 217 216 0 0 0 0 0 <NA>
## 218 217 0 0 0 0 0 <NA>
## 219 218 0 0 0 0 0 <NA>
## 220 219 0 0 0 0 0 <NA>
## 221 220 0 0 0 0 0 <NA>
## 222 221 0 0 0 0 0 <NA>
## 223 222 0 0 0 0 0 <NA>
## 224 223 0 0 0 0 0 <NA>
## 225 224 0 0 0 0 0 <NA>
## 226 225 0 0 0 0 0 <NA>
## 227 226 0 0 0 0 0 <NA>
## 228 227 0 0 0 0 0 <NA>
## 229 228 0 0 0 0 0 <NA>
## 230 229 0 0 0 0 0 <NA>
## 231 230 0 0 0 0 0 <NA>
## 232 231 0 0 0 0 0 <NA>
## 233 232 0 0 0 0 0 <NA>
## 234 233 0 0 0 0 0 <NA>
## 235 234 0 0 0 0 0 <NA>
## 236 235 0 0 0 0 0 <NA>
## 237 236 0 0 0 0 0 <NA>
## 238 237 0 0 0 0 0 <NA>
## 239 238 0 0 0 0 0 <NA>
## 240 239 0 0 0 0 0 <NA>
## 241 240 0 0 0 0 0 <NA>
## 242 241 0 0 0 0 0 <NA>
## 243 242 0 0 0 0 0 <NA>
## 244 243 0 0 0 0 0 <NA>
## 245 244 0 0 0 0 0 <NA>
## 246 245 0 0 0 0 0 <NA>
## 247 246 0 0 0 0 0 <NA>
## 248 247 0 0 0 0 0 <NA>
## 249 248 0 0 0 0 0 <NA>
## 250 249 0 0 0 0 0 <NA>
## 251 250 0 0 0 0 0 <NA>
## 252 251 0 0 0 0 0 <NA>
## 253 252 0 0 0 0 0 <NA>
## 254 253 0 0 0 0 0 <NA>
## 255 254 0 0 0 0 0 <NA>
## 256 255 0 0 0 0 0 <NA>unique(vat$Class) #classes for vat = various landscape types## [1] Background <NA>
## [3] Beach/Dune/Mud Open Water
## [5] Low Intensity Urban High Intensity Urban
## [7] Clearcut/Sparse Quarries/Strip Mines/Rock Outcrop
## [9] Deciduous Forest Evergreen Forest
## [11] Mixed Forest Row Crop/Pasture
## [13] Forested Wetland Non-Forested Salt/Brackish Wetland
## [15] Non-Forested Freshwater Wetland
## 14 Levels: Background Beach/Dune/Mud Clearcut/Sparse ... Row Crop/PastureWe can also see the number of cells in that class and the the refectances of the Landsat data at red, green, and blue spectral ranges. This is called a value attribute table (VAT). Notice that it runs from 1-256 (8 bits), but the size of the table does not correspond to the spatial extent, resolution, or the number of classes (18). This is quite different from a vector attribute table which includes a row for each feature.
Landsat Data = Data recorded by sensors of reflected and emitted energy.
Remember that we did a spatial join to characterize the landcover surrounding the detention basins? Let’s join this VAT to that shapefile so that the landcover values are more readily interpretable.
library(raster)
# Load detention basin points
dbasin_pts_utm<-shapefile("../IDEAS_Spatial_Ecology_Workshop/dbasin_pts_utm.shp")
# Create new field of the vat table named "lc" to join to the field of the same name in dbasin_pts_utm
vat$lc<-vat$Value
# Get only landcover value and description fields
vat2<-as.data.frame(vat[,7:8])
# Join tables
dbasin_pts_utm<-merge(dbasin_pts_utm,vat2,all.x=T) #joining by common names
# View results
head(dbasin_pts_utm)## lc X Y OBJECTI BasnTyp BasnSrf DatCrtd CnstrcT DamHght
## 10 22 -84.35676 33.88007 2801 Wet Pond Surface 1996 Riser 8
## 11 22 -84.35658 33.88119 2802 Dry Pond Surface 2007 Weir Wall 5
## 12 22 -84.35635 33.88207 2803 Dry Pond Surface 1989 Weir Wall 6
## 13 22 -84.35419 33.88468 2804 Dry Pond Surface 2001 Riser 10
## 1236 42 -84.35284 33.88555 2805 Dry Pond Surface 1989 Weir Wall 3
## 15 22 -84.35674 33.88303 2806 Wet Pond Surface 1999 Riser 5
## Age d_.HU_1 Class
## 10 23 Nancy Creek Low Intensity Urban
## 11 12 Nancy Creek Low Intensity Urban
## 12 30 Nancy Creek Low Intensity Urban
## 13 18 Nancy Creek Low Intensity Urban
## 1236 30 Nancy Creek Evergreen Forest
## 15 20 Nancy Creek Low Intensity UrbanIMPORTANT: For joins to the attributes of vector data, do NOT call the table as dbasin_pts_utm@data!! My guess is this overwirites the data stored in the original file? Or it breaks the entire item?
I just given a brief introduction to raster and vector data. In the next two modules, we will explore manipulating both in more detail.